Self-organized criticality - student project post-mortem

Comments

We just got finished with our student team project, which you can find here, and I thought I'd do a little post mortem on it.

It's a neat little project that implements various simulations of self organized criticality on 2D grids. What is self organized criticality, you might ask? Dunno, I can't tell you.

All right, I do know a little. Imagine the ising system I've written about before. In the version without an external magnetic field, it has a single important parameter that we can set - the temperature. If you sweep through the values of temperature, you can find a single point where the behavior of the system changes qualitatively - order wins over disorder at temperatures below roughly 2.72 in reasonable [set everything to 1] units. Near that value - at criticality - you get large scale behavior, huge fluctuations, exponential slowdowns, etc.

Self organized criticality, as I currently understand it, is basically that, except that as you run your simulation, you realize that it displays criticality for a wide range of parameters (here - temperature). To quote Wikipedia (emphasis mine):

the complexity observed emerged in a robust manner that did not depend on finely tuned details of the system: variable parameters in the model could be changed widely without affecting the emergence of critical behavior: hence, self-organized criticality.

Here's an example, a simulation of forest fire (yellow being fire, green being trees, and purple being ash, from which trees can regrow):

In [30]:
import SOC
model = SOC.Forest(L=100, f=0.0001)  # f being probability of lightning strike
model.run(500, wait_for_n_iters=500)
model.animate_states(notebook=True)

Waiting for wait_for_n_iters=500 iterations before collecting data. This should let the system thermalize.










    































What we did look for from a practical standpoint in our simulations was power law scaling of the number of iterations for each avalanche in a system, of avalanche size... and we did find that!


Here's another model we implemented, called the Manna model:













In [19]:



    


manna = SOC.Manna(L=30)
manna.run(500, wait_for_n_iters = 500)
manna.animate_states(notebook=True)



    















    






Waiting for wait_for_n_iters=500 iterations before collecting data. This should let the system thermalize.










    






 
 














    















    































This deserves a bit of a longer runtime:













In [23]:



    


manna_long = SOC.Manna(L=30, save_every = 100)
manna_long.run(50000, wait_for_n_iters = 20000)



    















    






Waiting for wait_for_n_iters=20000 iterations before collecting data. This should let the system thermalize.










    






 
 














    
























In [26]:



    


manna_long.get_exponent(low = 10, high=100)



    















    

















    






y = 23317.107 exp(-1.3904 x)










    
Out[26]:








{'exponent': -1.3903903386497791, 'intercept': 4.367674673827714}

























You can find the project on GitHub here.


Apparently, SOC has applications basically everywhere, but I haven't dug into those yet. There's are two articles on my to-read list that relate SOC to ELMs in fusion reactors and turbulent transport. I might do writeups on those in a later post.








 




    

    
Update on the year 2019

        

            

            

            
            
            

                

        
    Comments


        



    

    


2019 has been a weird year.


I don't feel like it's been all too good for my mental health, as I've definitely experienced burnout, mayhaps even depression. Suffice to say, the lows were low. I feel a bit like I procrastinated all there was to procrastinate.


It wasn't all bad, though - I'd say it was right up there in my top 25.


In trying to get some more work-life balance, I did get back to playing StarCraft - I'd cut it out of my life previously, along with all gaming - but I figured I do need some of that competitive drive in my life, and I still think there's a lot it can still teach me. You wouldn't drop chess from your life if you were passionate about it, right? So far, it's been treating me all right.


At the start of December, I started a research software job at the Institute for Plasma Physics and Laser Microfusion here in Warsaw, Poland! This is also where I'll be writing and developing my master's thesis, related to Bayesian inference for plasma diagnostics on Wendelstein 7-x, using Python and more specifically PyMC3. A step forwards for modern, open, reproducible and maintainable science software! At least, I hope so. :)


I was planning on putting some photos, stats etc here, but I haven't found the time to get to do those yet (remember all that procrastination? Yeah, it's coming back to bite me in the predictive posterior), and Matthew Rocklin says I should write shorter blog posts more frequently, anyway, and that sounds like good advice.


Thus, to  get this out the door more quickly... onwards to 2020!



    

    
Simple Binder usage with Sphinx-Gallery through Jupytext

        

            

            

            
            
            

                

        
    Comments


        



    

    


It's been a busy week for PlasmaPy. I recently found out about Binder support in sphinx-gallery. The latter is a package that we use to
turn python scripts with comments into Sphinx pages and Jupyter Notebooks. I figured adding that could be a nice fit for our existing example gallery .


However, I quickly realized that the system in place is a bit unwieldy. Binder takes a link to an existing GitHub repository and executes .ipynb notebooks located there online. However, with sphinx-gallery, we don't have those notebooks in the repository - we have .py files with comments. The currently recommended way of setting this up with sphinx-gallery is keeping your built documentation in another repository and hosting it via something along the lines of GitHub Pages rather than ReadTheDocs, which we are currently using.


I added the results of this investigation to sphinx-gallery's docs, but I didn't want to switch away from RTD, so I figured I'd go ahead and find another way. I think I've got something that works well enough now!


Trigger warning: later on during this post, there may be monkeypatching of sphinx_gallery internals. Beware.


Read more…



    

    
On the recent "On the Boris solver in Particle-in-cell simulations" paper

        

            

            

            
            
            

                

        
    Comments


        



    

    










I recently came across a pretty cool paper by Zenitani and Umeda named "On the Boris solver in particle-in-cell simulation". There are many splendid descriptions of the Boris solver on the Internet, so while I would rather not duplicate them, here's a brief overview. In PIC simulations, the Boris solver (or pusher) is the usual algorithm of choice for moving and accelerating particles in given electric and magnetic fields.


You may wonder, since the equations of motion are ordinary differential equations, what's wrong with using the usual Runge-Kutta 4 solver? As it turns out, that one has a pretty major flaw. It has great accuracy for short term calculations, but over time your particle's motion will lose energy. This is a deal breaker for periodic motion, and simulations of, for example, plasma waves need to conserve that energy to provide accurate results.


Boris came up with his solver in the 1950's, and in a single sentence: the algorithm splits the acceleration via electric field into two parts and sticks a rotation about the magnetic field between them. This turns out to conserve energy and will probably come up again on this blog as I read more about symplexicity. 



Read more…







    

    
Particle in Cell methods

        

            

            

            
            
            

                

        
    Comments


        



    

    


I think it might finally be about time to do some plasma physics
discussion on this blog, stay true to the name and so on…


Basically the only actual “scientific” work I have actually done with
plasmas up until now is writing a PIC simulation, PIC standing for
Particle-in-Cell. I thought I would take this opportunity to explain in
my own words what the concept is - I think it’s a clever one.


Read more…



    

    
Parsing and plotting LaTeX expressions with SymPy

        

            

            

            
            
            

                

        
    Comments


        



    

    












Today let's look into some pretty neat SymPy functionality. I was in a fluid dynamics lecture, practicing taking notes with LaTeX on the go and stumbled upon this monstrosity:

$$ \Delta(k) = \frac{\rho_1-\rho_2}{\rho_1 + \rho_2} gk + \frac{\gamma k^3}{\rho_1 + \rho_2} - \frac{\rho_1 \rho_2}{(\rho_1 + \rho_2)^2} U^2 k^2 $$
(bonus points for whoever recognizes this!)


We were supposed to draw this for a few example sets of values. All right! I opened up pinta and scribbled a few squiggly lines with my small touchpad, following the blackboard drawings. It looked darn ugly, but that got me thinking. SymPy has parsers, right? Can't I just parse that LaTeX equation into Python and make that plot pretty with matplotlib?


Well, as it turns out, sure...













In [17]:



    


the_plot.show()



    















    
































But it takes some tinkering.




Read more…